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ABSTRACT 



Of 



Context. The observations of rapidly rotating stars are increasingly detailed and precise thanks to interferometry and asteroseismology; two- 
dimensional models taking into account the hydrodynamics of these stars are very much needed. 
Aims. A model for studying the dynamics of baroclinic stellar envelope is presented. 
^ ■ Methods. This models treats the stellar fluid at the Boussinesq approximation and assumes that it is contained in a rigid spherical domain. The 
\Q temperature field along with the rotation of the system generate the baroclinic flow. 

f*^) . Results. We manage to give an analytical solution to the asymptotic problem at small Ekman and Prandtl numbers. We show that, provided the 
Brunt- Vaisala frequency profile is smooth enough, differential rotation of a stably stratified envelope takes the form a fast rotating pole and a 
slow equator while it is the opposite in a convective envelope. We also show that at low Prandtl numbers and without /^-barriers, the jump in 
viscosity at the core-envelope boundary generates a shear layer staying along the tangential cylinder of the core. Its role in mixing processes is 
O ' discussed. 

f~i Conclusions. Such a model provides an interesting tool for investigating the fluid dynamics of rotating stars in particular for the study of the 
C/3 various instabilities affecting baroclinic flows or, even more, of a dynamo effect. 

a ■ 
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1- Introduction 

Thanks to the development of new observational techniques, like interferometry or high precision photometry, rapidly rotating 
stars are more and more focusing the interest of the scientific community. The best exampl e is certainly the nearby star Altair 
whose shape, rotation, inclination of axis have been determined by in terferometry Jde Souza et all I20051 IPeterson et allEoOfjl) 
and which has also been identified as an oscillating 5-Scuti star (see iBuzasi et all I2005I) . Modelling such stars is therefore a 
challenge which needs to be taken up in order to extract the best science from these observations. 

The role of rotation in the physics of stars appears at various levels. First is the determination of the shape and thermal 
structure which control the way we "see" the star and at this stage the differential rotation is certainly important. Then, if we 
are concerned by the eige nspectrum, in addition to the shape (which influences the frequency and the visibility of the modes, 
seelLi gnieres et 3E005)), we also need to know the distribution of elements able to excite eigenmodes through K-mechanism. 
Finally, as it has long been known (e.g. ISpiegel & Zahnl Il970l) . rotation may drive many hydrodynamical instabilities in the 
stably stratified radiative zones of rotat ing s t ars w hich inevitably lead to some small-scale turbulence. Much work has already 
been done i n this direction following IZahnl dl992l) . in particular to understand the surface abundances in stars across the HR 
diagram (see Maeder & Mevnett l2000). 

The case of rotation is therefore fundamental in stellar physics and needs to be well understood. Presently, the difficulties 
concentrate on the hydrodynamical effects generated by the associated inertial accelerations (Coriolis and centrifugal ones). 
Indeed, the analysis of these flows is quite demanding as it requires two-dimensional models. 

The present paper focuses on one side of the general problem of the dynamics of rotating stars: namely, we want to specify 
the shape of the baroclinic flow inside the radiative envelope of a rapidly rotating star. In such stars, it is well known that the 
Eddington-Sweet time scale is short enough for a steady state to be reached. But this steady flow, in which viscous effects need 
to be taken into account, has never been computed in its full two-dimensionality . 

Previous work on this subject is quite scarce. First attempts include those of lTassoul & Tassoull (Il982lll983alblll984i ri986) 
who, for analytical tractability, restricted the velocity field to its first Legendre polynomial component and therefore give an 
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almost one dimensional description of the flow. Still in the quasi-one-dimensional approach, the work of lZahnl ll 19921) presents 
a self -consistent description based o n the assumption of a strong horizontal turbulent viscosity compared to the vertical one. 
Later, lUrvu & Eriguchillll994 [l995) touched upon the question of modeling rotating baroclinic stars in two dimensions but their 
neglect of viscosity removed any meridian circulation and left differential rotation incomplete. Including viscosity, Garaud ( 2002) 
investigated the dynamics of the radiative core of the sun stressed by the differential rotation of the convective zone; in this latter 
case the flow is essentially driven by the boundary conditions. This is certainly the most advanced work on stellar baroclinic 
flows but baroclinicity is not the main driving. 

Although t he full solutions were still not known, the stability of some generic situation of baroclinic shear flows has been 
investigated bv lKnobloch & Spruitl(H98lll983l) . lRnoblocrJ (1 1 985t) . ISchneideil ( 1 1 99dl) and are reviewed in lZahnl dl993l) . 

As one may easily imagine, computing the baroclinic flow in a centrifugally distorted envelope with realistic density and 
temperature profiles is not an easy task. Therefore, the step taken here concentrates on a simpler version of such a flow, namely 
a "Boussinesq version in spherical geometry", which, as we shall see, still contains a good deal of the actual physics of the 
problem. This model has two purposes: first to reveal the main dynamical features of such flows and eventually find some robust 
property that may persist in actual stars; second to devise a model in which hydrodynamics can be investigated more easily and 
which can be used as a template for the construction of more elaborated models. 

The next section of the paper describes our approach, especially the simplifications that we use or the physics we keep. 
In section 3, we discuss the asymptotic properties of this model in the limit of small Ekman numbers, appropriate to stellar 
applications. Then, after describing the numerical method, we investigate some aspects of the general case. Conclusions and 
outlooks follow. 



2. The Model 

2.1. Description 

We consider a system in which a self-gravitating fluid of constant density is enclosed in an undeformable sphere of radius R. The 
gravitational field is thus simply g = -gr where the radial coordinate is scaled with R (i.e. < r < 1) and g is the surface gravity. 
The thermal and mechanical equilibrium of this fluid is governed by: 



( -VP eq +p eq s = 

v ■ Orvr eq ) + q = o (i) 

[ Peq = P0(1 - ff(Teq ~ T Q )) 

where a is the dilation coefficient,^ the thermal conductivity and Q heat sinks (since we study a stably stratified situation). The 
essential characteristic of this equilibrium is the Brunt-Vaisala frequency profile 

N 2 (r) = a— ±g(r) (2) 
dr 

We now let this system rotating at the angular velocity Q. around the z-axis. In the rotating frame, because of the combination 
of rotation and stratification, a steady baroclinic flow appears. It is a solution of 



p(2Q A v + v ■ Vv) = -VP + p(g + n 2 se s ) + pAv 

p Cp v ■ vr = v • Orvr) + q 



where p is the dynamical shear viscosity, c p the specific heat capacity at constant pressure, s the radial cylindrical coordinate and 
e s the associated unit vector. We decompose the thermodynamical quantities into their equilibrium and fluctuating parts, namely 

P = P eq + 6P, p = Peq + dp, T = P e q + ST 

After subtracting Q, the momentum and heat equations read: 



p(2Q A v + v ■ Vv) = -V6P + 5p(g + Q 2 se s ) + p eq Q 2 se s + /xAv 
pc P (v ■ VP eq +V-V6T = V ■ CrV5P) 

These equations are further simplified by using the Boussinesq approximation which yields: 
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2fi A v + v ■ Vv = -V6P - a6T((g + Qrse s ) - a6T eq Q, 2 se s + vAv 

v ■ vr eq + v ■ vst = kMT 

At this stage we need to point out that we removed the barotropic contribution of the centrifugal acceleration, poQ. 2 se s , and 
that fluctuations of density are retained whenever they are multiplied by an acceleration (g or Q. 2 se s ) as should be done when 
using the Boussinesq approximation. We note that this approximation implies the Cowling approximation (i.e. the neglect 
of variations of self-gravity). We also simplified the system by assuming a constant heat conductivity. Taking the curl of the 
momentum equation we derive the equation of vorticity: 

V x (2£2 A v + v • Vv + a5T{{g + £l 2 se s ) - vAv) = -sN 2 (r) sin 6»cos Oe^ (3) 

where we used the Brunt-Vaisala' frequency profile (f2Jl and introduced s = Q?R/g, the ratio of centrifugal acceleration to gravity. 
These equations are of course completed by the equation of mass conservation which reads 

Vv = 

at the Boussinesq approximation. 

The foregoing equations show that our problem is that of a forced flow driven by the baroclinic torque -sN 2 (r) sin 9 cos Oe^. 

We may wonder at this stage how such a system compares with a real star and especially its baroclini city. In a star, the 
isothermal and isentropic surfaces are more spherical than the equipotentials or isobars (e.g. iBusseL Il982h : since the stellar 
envelope is stably stratified, entropy increases outwards which implies that, on an equipotential, entropy also increases from pole 
to equator. In our Boussinesq model, entropy (actually potential temperature) is represented by temperature; our equipotentials 
being oblate ellipsoids and the temperature gradient being positive outwar ds, we see th at there too, entropy increases from pole 
to equator of an equipotential surface. We may also note that as shown bv lZahnHl974l) . the baroclinic torque is proportional to 
the latitudinal gradient of entropy on an equipotential. 



2.2. Scaled equations 

Gathering vorticity, energy and continuity equations, we need to solve the following system: 
V x (2S2 A v + v ■ Vv + a5T(g + Q. 2 se s ) - vAv) = -sN 2 (r) sin 6 cos Oe^ 

v ■ Vr eq + v ■ V6T = kAST (4) 
Vv = 

since this is a forced problem, we wish to have a forcing of order unity as well as the solution. These considerations lead us to 
the following scaling of the velocity field and temperature perturbations: 

Q.N 2 R 2 , aT,g 

v = u, ST = sT t 8 with N = — — 

2g R 

where N is the scale of the Brunt-Vaisala frequency. This scaling takes into account the fact that the baroclinic flow and the 
associated temperature perturbations vanish when either rotation or Brunt-Vaisala frequency vanishes. We thus find the system 
of dimensionless dependent variables: 

( V x (c, A u + Ro u ■ Vm - (re r - sse s ) - EAu) = -tij(r) sin 9 cos Oe^ 

| (n 2 T /r)u r + su-V9 = E T A6 (5) 
[ Vh = 

where we introduced the numbers: 

2Q\ 2 N 2 R 



E T = T I ~TT I s RO 



2DJi 2 ' 2QJi 2 \N] ' Ag 

E is the Ekman number which measures the viscosity, Ej measures heat diffusion and Ro is the Rossby numb er. In addition t o 
these numbers we will need the Froude number, Fr, the Prandtl number V and the /l-parameter introduced bv lGaraud (E)02): 
these are respectively: 

V CINR _ v E 

Fr = = , r = -, A = — = P — t 

NR 2g k E T 4Q 2 

Finally, rij(r) is the scaled Brunt-Vaisala frequency. 
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Fig. 1. (a) Profiles of Brunt-Vaisala' frequency issued from the ID evolution code CESAM f or a 3 M n star at ( solid), 20 Myrs 
(dotted line), 284 Myrs (dashed line); these models include microscopic diffusion of Mic haud & Proffittl il993b . (b) Adopted 
profiles for a star with (dashed line) and without (solid line) yu-barrier. 



2.3. Boundary conditions 

This systems needs to be completed by boundary conditions. We shall assume the regularity of the solutions at the sphere's centre 
and impose stress-free boundary conditions on the velocity field at the outer surface. Thus doing, the velocity field is determined 
up to a solid rotation; if u is a solution of © then u + Ae z X r is also a solution (A is an arbitrary constant). For actual stars, such a 
degeneracy is lifted by initial conditions and conservation of angular momentum. Here we lifted it by imposing that the solution 
h of Q has no total angular momentum, i.e. that 

I su^dV = 
J 00 

so that the total angular momentum of the star is in the background solid rotation measured by Q. (which therefore appears as the 
mean rotation rate). 

We further complete Q by also imposing zero temperature fluctuation on the outer surface. 

2.4. Discussion 

As can be seen, the problem is controlled by a large number of parameters namely 



- the ratio of the centrifugal acceleration s to surface gravity which is also the ellipticity of equipotentials, 

- the Rossby number Ro, 

- the diffusion coefficients, E, E T , 

- the profile of Brunt-Vaisala frequency rij(r). 

Moreover, two other parameters will be necessary to described molecular weight gradients, namely the Brunt-Vaisala fre- 
quency profile n 2 (r), and the related diffusion coefficient, while another one will characterize the viscosity jump at the core- 
envelope interface. We therefore need some guidance in this large parameter space. 

For this purpose and in order that this model enlighten us on real systems, we shall consider the case of a 3 M Q star with a 
one day rotation period. Thus we will use a radius R = 2 10 9 m, a typical Brunt-Vaisala frequency of N = 10~ 3 Hz and a surface 
gravity of g = 10 2 m/s 2 . Besides, the profile of the Brunt-Vaisala frequency needs to reflect the more realistic models. In fig.^ 
we plot such a profile for a 3M Q star at different stage of its evolution on the main sequence. Such profiles should not be taken at 
face value especially t he cont ribution of the //-gradients since only microscopic diffusion is included in this model (produced by 
the code CESAM, see Morel ( 1997)). In fact it is the aim of the present work to understand the mechanisms by which elements 
move in the radiative envelope. Therefore we shall rather consider the generic situation visualized in Fig.^. 
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Fig. 2. (a) Profiles of the Rossby number rN 2 (r)/4g(r) issued from the ID evolution code CESAM for a 3M Q star at (solid), 
20Myrs (dotted), 284Myrs (dashed). 



Because of our choice of scalings, amplitude of nonlinear terms in the momentum equation is independent of rotation and 
can be appreciated directly from a non-rotating stellar model. In Fig. [2] we display the values of the Rossby number as a function 
of radius for our typical 3M Q star. Indeed, the non-dimensional function 



Ro(r) = 



N 2 (r)r 
4*(r) 



gives an interesting approximation of the amplitude of these non-linear terms. As can be seen, this number is generally less than 
unity except near the surface layers where it reaches Ro ~ 5. 

Recalling that e < 1, we see that the other non-linear terms are also less than unity (in fact very much less than unity as 
shown below). As a first step in the investigation, we shall ignore them altogether so as to be able to examine the properties of 
the system with linear equations. 

Setting Ro = s — 0, we obtain the system: 



V x (e z A u — Or — EAu) = —riL sin 9 cos 9e ¥ 



(nl/r)u r = E T A9 



V-h = 



(6) 



3. The asymptotic analysis at E «: 1 

Before solving the full system (|6j, we shall first discuss the case of asymptotically small Ekman numbers which are met in stellar 
applications. 



3.1. The inviscid profile 

When viscosity is neglected the equations admit a particular solution called (in geophysics) the thermal wind; in our case it 
reads 



u° = (sf^dr + F(s))e v , 
9 = 



(7) 



where s — r sin 9 is the radial cylindrical coordinate and F(s) an arbitrary function describing a pure geostrophic solution. 
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This solution deserves some comments. First, as already pointed out bv lBusseHl98lUl982l) . the so-called thermal disequilib- 
rium induced by centrifugal acceleration does not imply meridional circulation. A slight differential rotation where the baroclinic 
torque is balanced by the Coriolis torque gives a steady solution. 

However, we also see that such a solution is largely under-determined since the function F needs to be specified; this de- 
generacy is lifted by viscosity. The qualitative importance of F(s) comes from the fact that it controls the latitudinal differential 
rotation. 



3.2. The Ekman boundary layer 

The foregoing thermal wind solution usually does not satisfy the viscous boundary conditions and thus a boundary layer correc- 
tion needs to be added. Let us review here its main properties. 

If m° is the inviscid, interior solution and u the bound ary layer correction such t hat u° + u satisfies both the flow equations 
and the boundary conditions, it is well-known (see Greenspan, 1969, Rieutord, 1997) that u verifies 

e r A u + iu = C(u°) exp -\fi\ cos #|) 

where f = (1 — r)/ Ve is the radially stretched coordinate and C a complex vector. In the boundary layer the flow is essentially 
tangential and we may write the full velocity field: 

u g + iUip = C exp cos#|) + 

where the constant C is such that the horizontal stress is zero; namely 

d lu g + illy \ 
dr\ r 

which yields 

c = c(0) = (i + oy|r(0) 

with 

= F(sin 6) - sin flF'(sin 9) - « 2 (1) sin 
VI cos 6»| 

where the prime indicates derivatives. From these results we find the meridian velocity in the boundary layer: 



Ug = A / -r(<9)(cos f + sin <f)e^ (9) 



with £, = ( VI cos0|/2. 

Such a flow however does not verify mass conservation. This is taken care of by the so-called Ekman pumping which yields 
the radial component: 

u r = EU{sm6)cos£;e-t (10) 
with 



U(s) = n 2 (l)(l + q(s)) + sF"(s) + q(s)(F' - F/s), with q(s) = 1 + 



,v 2 



2(1 - s 2 ) 



These expressions show us that near the surface there exists a latitudinal flow 0(yfE) which induces an 0(E) circulation in 
the bulk of the fluid; in turn, this circulation controls the geostrophic flow F(s)e v and thus removes the degeneracy of the thermal 
wind solution. The reader may have noticed that the boundary layer solution is singular at equator (9 = n/2); this is the classical 
equatorial singularity of Ekman layers whose thickness changes to 0(E 2 ^ 5 ) in a region 0(E 1 ^ 5 ) around the equator (see also 
lGreenspanLfl969l) . 
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3.3. The thermal wind 

Since the differential rotation is a major feature of the baroclinic flow, it is useful to push further its analysis and try to express 
F(s) as a function of the Brunt-Vaisala' frequency profile. 

For this purpose, we need to express the meridian circulation as a function of F. We thus write the ^-component of the 
momentum equation in cylindrical coordinate 

u s = E(A-l/s 1 )ul 

where A is the laplacian operator; this equation expresses the local balance between advection and diffusion of angular momen- 
tum. It yields 



u s = E(F" + F'/s - F/s 2 + sC(r)), 



with 



, .. , An 2 
C(r)= - +-r 



Mass conservation gives the z-component of the flow: 



u, = — 



1 dsu, 



o s ds 



dz = -EL 3 (F)z 



-Efi 
Jo s 



ds 2 C(r) J , 

a dz 

OS 



(ID 



(12) 



with L 3 (F) = [(sF')" - (F/s)'] js while in the integral r 2 = z' 2 + s 2 . 

Now, we put all the pieces together by demanding that u r — is satisfied at order 0(E) on the outer boundary. Hence, setting 
r — 1 , we have 

su s + zu z + u r — 

which gives the following differential equation for F: 

C ({s) I C'(r)\ 

sF" + F' -F/s + s 2 C(l)-{(s) 2 L 3 (F)-((s) J \2C(r) + s 2 —-^-\dz + n 2 (l)(l + q(s)) + sF"(s) + q(s)(F' - F/s) = (13) 

with £(s) = Vl - s 2 . This complicated equation is of the form -C(F) = b, namely a forced third order differential equation. We 
shall not try to find the general solution of it, but will examine the polynomial solution F = As + as 3 + bs 5 , formally valid at 
s <k 1, with the hope that since s < 1, it will give a rather good idea of the flow. 

We first note that, as expected, the solution is invariant to the addition of a solid rotation; thus A is arbitrary. The first 
interesting term is as 3 . Substitution into (II 3i give the following result: 



1 n 2 (r) 



dr, 



b = 



1 In 



48 \ r 



(1)+ l92 n(1) -24j U + 3r[[- 



dr 



and thus we expect the following differential rotation close to the z-axis 



C n 2 (r) , s 2 f 1 n 2 (r) 

J -r dr -2j — 



(r) dr + bs 4 



r ^ Jo ^ 

which is now an explicit function of the Brunt-Vaisala frequency profile. 



(14) 



3.4. The role of stratification 

In the foregoing solution we set 9 to zero and therefore did not bother about stratification. In order to understand its role, we need 
to evaluate the amplitude of 9, i.e. the scaled temperature fluctuation. 

Equation (|Sp) shows that since u r (i.e. meridian circulation) is 0(E), then 9 is 0(E/Ej) that is of the order of A — VN 2 . Thus, 
in the limit of small A's, our solution is certainly valid. Stellar fluids are characterized by small Prandtl numbers and usually in 
rapidly rotating stars one considers that A <k 1 which makes our neg lect of buoyancy comfortable. 

However, as was pointed out long ago by Spie gel & Zahnl Jl970l) it is most likely that some turbulence is driven by the shear 
of the baroclinic flow. The viscosity, change d to a turbu lent one, is then increased by some factor which therefore also increases 
the Prandtl number. From the discussion of Zahn ( 1993), it turns out that this turbulent Prandtl number may reach values close 
to unity; hence, even in a rapidly rotating star, the /l-parameter may be larger than unity (typically up to 10 2 ). 

The consequence of such a situation is that a thermal boundary layer appears. Its scale is of course 1/ Viand we now assume 
that A » 1. Outside this layer, quantities vary on an 0(1) scale; noting that u s ~ u r , using the ^-component of the momentum 
equation together with the heat equation, it follows that 

6 ~~ Alitn 
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Fig. 3. Differential rotation (right) and meridional circulation (left) for a Brunt-Vaisala frequency profile n(r) = r. 

Now, the ^-component of the vorticity equation (|6^) shows that in this region 

~ A' 1 , u r ~ EA , and 6~1 

Roughly speaking, outside the thermal layer the temperature fluctuation is of order unity and the velocity perturbation is reduced 
by a factor 0(A), while in the thermal layer, radial gradients are increased by a factor Va which gives: 

u v ~ A , u r ~ E, and 6 ~ 1 

Hence, all the dynamics described above for P« 1 becomes squeezed inside the thermal layer with a differential rotation reduced 
by a factor A. This is in fact very similar to the solar tachocline where turbulent diffusion inhibits the diffusion of momentum 
with the help of stratification. 

3.5. Discussion 

The foregoing solution potentially shows very interesting properties of the baroclinic flow; however, we may wonder if it is 
reliable mathematically considering the crudeness of the solution of Jl 31 . 

The accuracy of the as 3 solution of Jl 3i may only be estimated by direct comparison to numerical solutions. Anticipating on 
the next section we solve numerically l|6) in the limit of small Prandtl number. We used a very simple Brunt- Vaisala frequency 
profile, namely n 2 = r 2 , which is actually quite standard in the literature (see Chandrasekhar, 1961). With such a profile the 
estimate of the different components of the flow are quite easy to evaluate (a = — | and b - j^); as for the differential equation 
we find 

sa = \t 

neglecting b; this solution compares extremely well with the numerical solution (see Fig. [3} as soon as P < 0.1. It shows that the 
second term bs 5 does not much influence the shape of the solution. From the expression of b we see that this is true because, in 
this case, n 2 is a smooth function of r. In stellar situations, the coefficient (« 2 /r)'(l) might be very large if sharp gradients develop 
in the upper layers of the star. 

Reassured by these computations, equation MAI therefore shows a very interesting result: the Brunt-Vaisala frequency 
profile of a star is sufficiently smooth, the latitudinal differential rotation in a stably stratified envelope driven by the sole baro- 
clinicity is a fast rotating pole and a slower rotating equator. Moreover, the a-coefficient being an integral over the star, it is not 
sensitive to the detailed shape of the Brunt-Vaisala profile and therefore is a rather robust quantity. 

Conversely, in a convecting envelope, with viscous-like Reynolds stresses, baroclinic torques induce a fast rotating equator 
and slow poles, like the sun actually. Of course, in the solar case, this may just be a coincidence since baroclinicity in the 
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convective zone of the sun hardly comes from centrifugal effect but rather from the convection itself (see lBrun & Toomrell2002l) . 
Nevertheless, in rapidly rotating (i.e. young) solar type stars, our result shows that baroclinicity drives a solar-like differential 
rotation. 

Finally, the solution also gives the form of the meridian circulation through dl It and fl!2i . For instance the number of cells 
can be retrieved from the expression of u s . At equator, s — r and u s reads 

/ 1 d(Pn 2 ) f 1 n 2 

Us = E \p— r 4r 1 ~ ( ' r 

The number of zeros of this function, plus one, gives the number of cells according to the Brunt-Vaisala frequency profile. 

Besides, the expression of ue in the Ekman layer © shows that at the very surface, ug n 2 (l)s, which means that the flow is 

poleward. If the Brunt-Vaisala frequency profile is smooth enough this motion extends to low latitudes. However, because of the 

sin £ + cos £ dependence, the boundary layer flow reverses just below the surface at r = 1 — ^ yj^fgj. 

From the amplitude of the meridian flow, one can estimate the circulation time scale which is: 



2Q T. 



vise 



■ circ « r J-. 1 vise t-. 

NFr Ro 

Since Ro ~ 1, we see that this time scale is of the order of the viscous diffusion time scale 7Vi SC which is very large if only 
microscopic viscosity is diffusing momentum. 



4. The general case 

We now turn to a more complex model with which we wish to analyse the effects of a convective core and of yu-gradients. The 
inclusion of yu-gradients amounts to the addition of the equation of a concentration c and the modification of the buoyancy term; 
the linear system (|6} needs to be changed into: 

V x (e z A u - (6 - c)r - EAu) = -n 2 sin 9 cos 6e 9 
{n 2 T lr)u r = E T A6 

(15) 

-{n^lr)u r = E c Ac 
{ V • u = 

where we introduced n 2 = nj + n 2 and E c = {D C /2QR 2 ) x (2Q/W) 2 with D c being the diffusivity of element "c". This element is 
supposed to be heavier than the surrounding gas. 

Since solutions are much more complicated we turn towards numerical solutions and first describe the numerical technique. 

4.1. The numerical method 

To solve the system (I15> we first project the variables onto the spherical harmonic base (see iRieutordi 1 1 9 87l) 

« = ZZ u >" R ™ + v <« 5 "' + w - r ™ 9 = Z Z ^ Y "'' c = Z Z c >» Y "' 

e=o m=-t e=o m =-e c=o m =-e 

where 



R m t = Y™e r , = VYf, T" e ' = Vxfi 



with Y? being the usual normalized spherical harmonic function. Since the flow is divergenceless, v' m may be expressed as: 



e 

1 1 dr 2 ul 



1(1 + 1) r dr 
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Then, projecting the equations (the vorticity equation needs be projected only on R™ and T™), we get 

EA ( w ( m + A{()r^ + A{( + ly-^^ui^ = 

EA ( A ( {ru e m ) = Btf)/- 1 ||$r) + B{€ + l)^ 2 ^ 2 ^ 1 J + t{l + \)((f m - c c m ) - {n\ + tfyN 2 8 a 



E T A e 6 e m - n\lru l n = 
{E c A t ci+nllmi l = Q 

where N 2 = J±¥ and 



P \{2t- \){2£+ 1) ' y (2£ - l)(2£ + 1) 

This coupled system of differential equation needs to be completed by boundary conditions at the surface of the sphere (at 
the centre we only demand the regularity of the solutions). Concerning the velocity, we impose stress-free conditions, thus 

while on the temperature and concentration we impose respectively no fluctuation and no flux at r = 1, namely 

<-«. § = ° 

Further, we shall need interface conditions at the core-envelope boundary. These conditions express the continuity of the 
velocity, temperature, concentration, fluxes and stress fields. The continuity of first quantities simply translates as the continuity 
of 

u m> -^T> W m' U m' c m { - lb > 

while the continuity of fluxes and stresses must reflect the changes in the transport coefficients between the turbulent convective 
core and the radiative envelope. Continuity of the three components of stress demands the continuity of 



dr r 

dui / d 3 ui , d 2 uL „, dut 



rt-ZE-S or EW^ + 6 r-^ + 3(2-^+1)) 



<Hn\ 

dr J 



while continuity of heat flux and concentration flux impose the continuity of 

£ T d Jk and e/A 
dr dr 



4.2. The shape of the baroclinic flow 

We first investigate the baroclinic flow which is represented by both the differential rotation and the associated meridian circula- 
tion. We focus our attention on the influence of the main uncertainties of the problem, namely the sensitivity to viscous transport 
and the height of the ^-barrier. 
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Nr=130 L=120 E = 1.0xl0-« Jfel.OxlO-* N 2 m ,„=1.00 v = 1.00 CL=ft Nr=130 L=120 E=1.0xl(H *>=1.0xl0-« N 2 m „„=1.00 i/ = 1.00 CL=ft 



Fig. 4. Left: the meridional circulation of the baroclinic flow generated by the Brunt-Vaisala frequency profile in figure^) (solid 
curve) for a fluid with constant viscosity. Right: The associated differential rotation showing the fast rotating pole and slow 
equator (solid lines are for positive contours, dotted lines for negative ones). 



4.2.1 . With no jump in viscosity 

We first consider a simple configuration where we set a Brunt-Vaisala frequency profile with no //-barrier, no viscosity jump at 
core-envelope interface. We thus find a plain baroclinic flow and can check the agreement with our foregoing theoretical results. 
The meridian circulation is shown in Fig.0]together with the associated differential rotation. The plot of the angular velocity show 
the typical shape of the thermal wind solution given by @. This shape remains identical if the viscosity is decreased (i.e. with 
lower Ekman or Prandlt numbers). This flow is essentially azimuthal with faster rotating poles as predicted above. The meridian 
circulation is much weaker (see Fig. |5jl being essentially 0(E). As expected from the above boundary layer analysis, near the 
surface the latitudinal component increases very strongly by a factor E~ 112 (this is essentially an effect of mass conservation). A 
detailed view of the velocity profiles in the Ekman layer is given in Fig. [6] 

To have further information on the baroclinic flow we had a look to the case where the Prandtl number is of order unity. In that 
case the value of A is important and we set it to 100 according to our typical rotating star. The meridian circulation along with the 
radial profiles of the velocity field are shown in Fig0 Obviously, the circulation has much decreased. This is a direct consequence 
of the heat equation (|6j>) which imposes a reduction of radial velocity when heat diffusion decreases; as a consequence advection 
of momentum is less and therefore diffusion of momentum must reduce which is obtained by reducing the differential rotation as 
shown in FigEJ>. 



4.2.2. With a viscosity jump at the core-envelope boundary 

In our quest of more realistic models, we now consider the effect of the viscosity jump at the core-envelope boundary which 
is brought about by turbulent convection in the core. The convective core is thus considered as a much more viscous fluid with 
negligible stratification (Brunt-Vaisala frequency is set to zero). 

The interesting result is that the meridian circulation is strongly modified as shown in Fig [8] This figure shows that the jump 
in the mean mechanical properties of the fluids generates a shear layer parallel to the axis of rotation. Such a layer is a Stewartson 
layer well-known in the dynamics of rotating fluids. The dynamics of these layers is controlled by a delicate balance between 
viscous stress, pressure gradient and the Coriolis acceleration. As shown long ago bv IStewartsonl i 1 96(h . such layers are in fact 
nested layers whose width scales with v 1 ^ 3 , v 1/4 or v 2/7 . However, if the properties of such layers are well known in the simple 
case of incompressible fluids, they remain unexplored when the fluid is stably stratified like here. 
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Nr=230 L=120 E=1.0xl0- 6 j»=1.0xl0- 4 N 2 ,„„=1.00 v = 1.00 CL=ft 



Nr=230 L=120 E=1.0xl0- 6 j»=1.0xl0- 4 N 2 „,„ =1.00 i/ = 1.00 CL=ft 



1.0 0.£ 



0.96 0.97 0.98 0.99 



Fig. 5. Radial profile at 9 — lrd of the velocity compo- Fig. 6. Same as Fig [5] but with emphasis on the Ekman 
nents of the baroclinic flow. Note that the v r and vg have layer, 
been multiplied by a factor 0(7? _1 ' 2 )» 1. 




Fig. 7. Left: Same as Fig.0t but with P — 1 . Right: The associated velocity profiles. Note the squeezed circulation in the thermal 
layer which width is ~ 0.1, and the reduction of the azimuthal velocity in the bulk of the domain. 



The interesting point of this feature is that the scale of the meridian flow is much reduced in one direction; the scale controlling 
the shear layer pumping is 0(E 1 ^ 4 ) if we follow Stewartson theory. This means that the circulation time is reduced by a factor 
0(E 1 ^ 4 ). Indeed, mass conservation implies that 

dV ± _ av„ 

dxj_ dx/i 
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where V± is the pumping induced by the shear layer and which should be identified to the meridian flow (see figure |9j. Since 
d/dx ± ~ E~ l t 4 and d/dx// ~ 1, it turns out that V// ~ E~ x ^ V±; hence the advection time from core to surface is reduced by a 
factor 0(E 1 ^ 4 ) which is quite small as shown below. 

Now we may wonder whether this dynamical feature will resist to the building of /i-barriers during the evolution of the star. 
To answer this question, we computed a model with a yu gradient as shown in Fig. [Q> (dashed line); namely the height of the 
/i-barrier has been raised such that ~ Nj. As far as elemental diffusivity is concerned we assumed that it is of the same order 
of magnitude as the viscosity, i.e. the Schmidt number v/D c is of order unity. Figure [TO] shows the result: the Stewartson layer 
has disappeared. It is replaced by a meridional flow within the core. 



14 



M. Rieutord: The dynamics of the radiative envelope of rapidly rotating stars 




Fig. 11. Angular momentum distribution and its derivative. The dotted lines in dL z /d8 show regions where the barotropic ax- 
isymmetric instability may develop (solid lines are for positive contours, dotted lines for negative ones). 



4.3. Questions of stability 

After the computation of the global baroclin ic flow, we may wonder w hether it is stable or not. This question is not an easy 
one since many instabilities are possible (see iKnobloch & Spruiti fl982i) and the investigation would require a dedicated paper. 
However, we can discuss qualitatively the question. 

A first class of instabilities are the barotropic ones; they do not take advantage of the baroclinic state and are basically due to 
the shear. The first to be considered is the axisymmetric one; it is also known as the centrifugal instability and is controlled by the 
angular momentum profile. The Rayleigh criterion, taking into account the stable stratification, says that instability sets in when 

dL, „ 1 dL 7 „ 

— - < at r = Cst, <=> < 

os cos 6 88 

that is when the angular momentum decreases with the cylindrical radius. The constraint r = Cst avoids the stabilizing effect of 
the Brunt-Vaisala frequency profile. 

For our system the angular momentum of the flow is 

L z = Qi 2 + sV v = OR 2 (s 2 + 2Ro su^) 
so that the flow is unstable if 

r cos 9 + Ro W + l ) w L Y "' < 

e 

This formula shows that the Rossby number controls this instability. As discussed above, this parameter is in the range [0.1, 10] 
in our reference star. We find that for Ro < 3 the flow is stable and, as shown in Fig. equatorial regions at destabilized first 
when Ro increases. We shall not push further the analysis since at this point only a detailed analysis would give sensible results. 
We just conclude that in view of Fig.|2] we expect such instability only in equatorial surface layers. 

The other barotropic instability is the non-axisymmetric shear instability. Basically, in a stably stratified fluid it is controlled 
by the Richardson number, e.g. 

N 2 , 

Ri = — t = Fr~ 2 

{dvpl 'dry 

In our case this number appears to scale with the inverse square of the Froude number which we found of order unity. We conclude 
that most likely the large-scale flow is stable and that a rigorous stability analysis is needed to give precise answers in stellar 
conditions. Nevertheless, following lZahnl \ 1 99 3). this flow may still destabilize scales small enough that they are insensitive to the 
stable stratification because of the large heat diffusion coefficient. The instability of small scales might thus provide the rotating 
radiative zone with some turbulence and hence some turbulent viscosity. IZahnN 19931) shows that such turbulent viscosity should 
be within the range 

Re c v.v r =, 2 -^ — (^) 

where Re c is the critical Reynolds number for shear flow, typically around 1000, and Ri ( . the critical Richardson number for stably 
stratified shear flow, around 0.25. These bounds on the turbulent viscosity show that the viscosity is increased at least by a factor 
Re e and at maximum by a factor 0(fP~ l ). Interestingly enough, this raises the Ekman number from 10~ 18 up to ~ 10~ 15 , 10~ 12 
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which is still very small. Moreover, the Prandtl number is increased by a factor between Re f and V which means, in the latter 
case, an effective Prandtl number of order unity and a /l-number significantly larger than unity. In fact, the turbulence is likely 
anisotropic and the values of the effective Prandtl number depends whether one considers horizontal or vertical diffusion. An 
effective Prandtl number of unity is therefore an upper bound. 

Besides these classical instabilities (for fluid flows), w e have also to face barocli nic instabilities which are specific to this 
type of flows driven by the baroclinicity. As discussed by ISoruit & Knoblocr] dl984h . one should separate the true baroclinic 
instability which develops because some fluid parcel may move in a direction between entropy levels and isobars, and diffusive 
instabilities (Goldreich-Schubert-Fricke and ABCD) which take advantage of double diffusion. These instabilities, specific to 
baroclinic flows, are rather weak and one question is whether they would persist if the diffusivities are increased by some 
turbulence generated by the barotropic ones. More work is needed to characterize them, especially in a global approach. 



5. Discussion 

5.1. Stellar numbers 

To fix ideas on the foregoing results we now give some numbers using our test star. We recall that it is a 3 M rotating with 
a period of one day. On the ZAMS, its radius is 2.6 Rq = 2 10 9 m. Typical kinematic viscosity of the plasma is v = 10~ 3 m 2 /s 
yielding an Ekman number of 

£~21(T 18 



really very small. The corresponding viscous time is r v ; sc = R /v ~ 10 yrs. The circulation time is typically a viscous time: 

rj. _ R 1 -^visc 

* ciic — 



V r FlN U<p Ro 

since Ro ~ 1 . As we discussed above, various instabilities may trigger turbulence and hence raise the viscosity up to turbulent 
values. From the typical value of the microscopic Prandtl number, namely 1CT 6 , and that of critical Reynolds numbers, ~ 10 3 , 
we can expect and increase of Ekman number up to 2 1CT 15 — > 2 10~ 12 and hence a reduction of the circulation time down to 
10 11 — > 10 8 yrs (10 8 yrs is the Kelvin-Helmoltz time of our model). In the early phase of evolution, when /i-barriers are not strong 
enough, this time scale may be further reduced by the presence of the Stewartson layer. However, this layer is sensitive to the 
/l-number, which means sensitive to the Prandtl number. Thus, if turbulence is strong enough to raise the Prandtl number to order 
unity values, the Stewartson layer disappears. On the other side if it is mild enough and increases the viscosity by a factor 10 3 
only, the Stewartson layer has a thickness 0(2 10~ 4 ) which reduces the circulation time further by a factor 510 3 ; all in all, mild 
turbulence and Stewarson layer may reduce the circulation time by a factor 510 6 leaving a time scale for partial mixing around 
10 yrs. On the other hand if turbulence is strong, circulation outside the thermal layer, which reaches the core, is ~ P~ l faster 
thanks to the increased viscosity, but A. slower because of stratification. However, the A~ l factor only reduces the advection 
time, if turbulence is strong, elements may diffuse just like momentum and the time scale is 7\,j sc !P = 10 8 yrs. 

The point raised here is, rather than the numbers, the role played by the different hydrodynamical processes in controlling the 
transport of element. Interestingly enough, we see that these processes may control each other: a mild turbulence being helped by 
a strong Stewartson layer yields a mixing time scale not much different from the one derived in the strong turbulence case. This 
leaves the possibility that mixing may not be too sensitive to the hydrodynamical details. In particular, we do not include angular 
momentum loss through a stellar wind (see lZahnlll992l) : our result shows the possible mixing occuring in the weak wind limit. 

To end this section let us put some numbers on the differential rotation induced by baroclinicity. The scalings give the 
following relation: 

AD. 

— = 2Ro(5Q 

where SO. = u^/s. The numerical solutions show that <5O p0 i e - <5Q e q ~ 0.5 when P « 1 and decreases to 0.01=O(A~ l ) when 
<P~\. 

From these values we see that such a model predicts a rather important differential rotation when the Prandtl number is small 
and a strong reduction of it if some turbulence develops. 



5.2. Non-linear terms 

Discussing the results, we may wonder about the importance of the non-linear terms that we neglected. They are Ro (u ■ V«) 
and su ■ V#. Neglecting 0(E) terms, the momentum nonlinear term is just -Ro u^/se^, i.e. a centrifugal correction to that of the 
background rotation. We thus expect a slight modification of the baroclinicity but no qualitatively important change. 
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Concerning the heat advection, the weakness of meridional circulation 0(E), plus the 0(A) amplitude of 6, make this nonlinear 
term O(sAE), thus extremely small and really negligible. 

It hence turns out that nonlinear terms are either very small or not bringing new physical phenomena in the steady solution. 
Taking into account the initial assumption of the Boussinesq approximation, refining the solution with these terms does not make 
sense. 

5.3. Back to the Boussinesq approximation 

To conclude this discussion, we briefly come back on the use of the Boussinesq approximation. One may indeed wonder 
whether in a more realistic approach, taking into account density variations of the background, our results would persist. 
The Boussinesq approximation is indeed the consequence of two limits: the fluid velocity is small compared to the sound 
velocity and the density scale height is large compared to size of the container (see Rieutord. 19971) . The first constraint is 
easily met: Mach numbers are small. The second one is not because the density varies on a scale comparable to the size 
of the star. Thus, the simplification of the equation of mass conservation from V ■ pv = into V ■ v = is certainly abusive 
and quantitatively our model should be taken with care. On the qualitative side, however, we note that density variations 
cannot modify the hydrodynamical features like the sign of the differential rotation, the appearance of Stewartson layer, 
etc. One may actually note that the use of the momentum pv instead of v would solve this problem, but then be inconsistent 
with the neglect of the effects of centrifugal distorsion. Thus, even if at first sight the use of the Boussinesq approximation 
seems exagerated for a star, it provides an interesting view of its hydrodynamics: it takes into account all the force fields 
which control the flows and gives a physically self-consistent model. 

5.4. Conclusions 

To end this paper we would like to stress the most interesting points of this model, namely: 

- The determination of the differential rotation in radiative envelope as a function of the Brunt-Vaisala frequency profile and 
its sensitivity to the amplitude of a turbulent viscosity, 

- The demonstration that, when the Brunt-Vaisala frequency varies on a scale of the order of the radius, baroclinicity generates 
a fast rotating pole and slow equator and the opposite in the case the envelope is convective, 

- The determination of the meridian circulation, its shape and amplitude, its sensitivity to Prandtl number and Brunt-Vaisala 
frequency profile, 

- The appearance of Stewartson layers induced by the jump in the mechanical properties of the fluid at the core-envelope 
boundary and the sensitivity of this layer to the /i-gradients or the Prandtl number, 

- The richness of the dynamics of these envelope and the compensation, as far as mixing is concerned, which may result from 
large-scale flows (Stewartson layer) and small-scale turbulence. 

Naturally, such a model cannot pretend to describe actual stars at face value. It should rather be viewed as a laboratory 
experiment aimed at studying some aspects of the stellar dynamics. Its (relative) simplicity indeed authorizes the study of hydro- 
dynamical instabilities at global scale through either linear analysis or full numerical sim ulations. For instance, one may think to 
use it to study dynamo effect in radiative envelope following lBraithwaite & Spruitl J2005 ) and to study the effects of an additional 
driving like the angular momentum loss generated by a wind. 

The next step is to compute the same flow in the real, spheroidal, geometry of rapidly rotating stars. Such realistic models 
will determine the aspects of the dynamics of stellar envelopes which are robust and can be studied with our Boussinesq model. 
They may also suggest some intermediate (artificial) models where a baroclinic torque is just plugged in a spherically symmetric 
stellar model for studying the dynamical aspects which are not sensitive to the true spheroidal shape. 
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